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ABSTRACT 

It was predicted more than 40 years ago that the cores of the coolest white dwarf stars should eventually crystal- 
lize. This effect is one of the largest sources of uncertainty in white dwarf cooling models, which are now rou- 
tinely used to estimate the ages of stellar populations in both the Galactic disk and the halo. We are attempting 
to minimize this source of uncertainty by calibrating the models, using observations of pulsating white dwarfs. 
In a typical mass white dwarf model, crystallization does not begin until the surface temperature reaches 6000- 
8000 K. In more massive white dwarf models the effect begins at higher surface temperatures, where pulsations 
are observed in the ZZ Ceti (DAV) stars. We use the observed pulsation periods of BPM 37093, the most mas- 
sive DAV white dwarf presently known, to probe the interior and determine the size of the crystallized core 
empirically. Our initial exploration of the models strongly suggests the presence of a solid core containing 
about 90% of the stellar mass, which is consistent with our theoretical expectations. 

Subject headings: stars: evolution — stars: individual (BPM 37093) — stars: interiors — stars: oscillations — 
white dwarfs 



1. MOTIVATION 



Mor e than four deca des have passed sin ce lAbrikosovl 
(1960), Kirzhnitz (1960) and Salp eteil (119611) predicted that 
the cores of white dwarf stars should crystallize as they cool 
down over time. There has never been a direct empirical 
test of this theory. The discovery of pulsations in the mas- 
sive hydrogen-atmo sphere (DA) white dwarf BPM 37093 
( Kana an et alj fl992) provided the first opportunity to search 
for the observational signature of crystallization in an indi- 
vid ual star. Theoretical calculatio ns bv IWinget et alJ ( 1997) 
and Montgomery & Winget ( 1999) suggested that the core of 
this star might be up to 90% crystallized, depending on its 
mass and internal composition. 

In addition to providing the first test of the theory of crys- 
tallization in a dense stellar plasma, knowing whether and 
to what degree this star is crystallized has implications for a 
more fundamental question. Recent Hubble Space Telescope 
observa tions of t he faintest white dwarfs in the globular clus- 
ter M4 ( lHansen et al.l2002l) have led to a resurgence of inter- 
est in using these stars to provide independent constraints on 
the ages of stellar populations. The largest potential sources 
of error in this method arise from uncertainties about the com- 
position and structure of white dwarf interiors. Fortunately, 
all of the major uncertainties can be minimized through de- 
tailed observation and modeling of pulsating white dwarfs, 
providing a crucial method to probe the stellar interiors and 
calibrate the cooling models. 

The crystallization process leads to one of the largest 
sources of unce rtainty in the ages of cool white dwarfs 
dSegretain et al.|1 994l). When a typical mass w hite dwarf star 
(0.6 Af ? . iNaniwotz ki. Green. & Safferlll999i) cools down to 
T e ff ~ 6000-8000 K (depending on the core composition), the 
high-density core will undergo a phase transition from liq- 
uid to solid. An associated latent heat of crystallization will 
be released, providing a new source of thermal energy that 
introduces a delay in the gradual cooling of the star (for a re- 
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cent review, see Han sen & Liebertl2003l) . In mixed C/O cores, 
phase separation of the ions during crystallization can provide 
an additional source of energy, delaying the cooling even fur- 
ther. 

Early attempts to determine the crystallized mass fraction 
in BP M 37093 were plagued by difficulties with uniqueness 
(Montgomery & Winget 1999). Recent improvements in our 
ability to match the observed pulsation periods in white dwarf 
stars with theoretical models have been driven by the devel- 
opment of an optimization method based on a parallel genetic 
algorithm (Metcalfe & Charbonneau 2003). This method al- 
lows the objective global exploration of the defining param- 
eters, which is essential to ensure a unique solution. In this 
Letter, we present the initial application of this method to fit 
the observed pulsation periods of BPM 37093 with asteroseis- 
mological models. 

2. OBSERVATIONS 

BPM 37093 has been the target of two multi-site ob- 
serving campaign s of the Whole Earth Telescope (WET; 
iNather et al.ll990l) . Preliminary results f rom these campaigns 
were published by Kanaan et al. (2000). The 1998 observa- 
tions (XCOV16) revealed a set of regularly spaced pulsation 
frequencies in the range 1500-2000 pHz, which iNittal (120001) 
identified as six different radial overtones (k), possibly hav- 
ing the same spherical degree (£ = 2). The 1999 observations 
(XCOV 17) revealed a total of four independent modes, includ- 
ing two new modes and two which had been seen in the pre- 
vious campaign. We obtained new single-site observations 
of BPM 37093 from the Magellan 6.5-m telescope on three 
nights in February 2003. These data showed evidence of five 
independent modes, all of which had been detected in the two 
previous multi-site campaigns. 

The Fourier Transforms (FTs) of the two WET campaigns 
and our new single-site observations are shown in Fig. \l\ 
Taken toget her, these data support the hypothesis — originally 
proposed by Kleinman et al. (1998) — that each set of observa- 
tions can provide a subset of the full spectrum of modes that 
are actually excited in the star. For reasons that are not well 
understood, the amplitudes of many of the pulsation modes 
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FIG. 1. — Fourier Transforms of the observed light curve of BPM 37093 
from multi-site Whole Earth Telescope campaigns in (a) 1998 and (b) 1999, 
along with (c) single-site observations from Magellan in 2003. Each panel 
also contains the corresponding window function — the pattern of alias peaks 
present in the Fourier Transform for each real frequency. The amplitudes of 
various pulsation modes change over time, but their frequencies are stable. 



observed in cool DAV stars are highly variable, even though 
the frequencies appear to be relatively stable. For example, 
the three FTs shown in Fig. ^ a H contain the two closely- 
spaced modes near 1575 /iHz. However, the mode near 1720 
/iHz in Fig.^ was not detected the following year in Fig.^3, 
but returned with the same frequency in Fig.^;. Likewise, no 
mode is detected near 1875 /iHz in Fig.^, while it is clearly 
present in both Figs.^i and^. 

We adopt the full set of 8 independent modes proposed by 
iNittal d2000l) for model-fitting. We as sume that each s ingle 
mode has an azimuthal order m = (see;Metcalfe 2003a]). For 
the three modes consisting of two closely-spaced frequencies, 
we use the average of the two. The adopted list of observed 
periods is shown in Table [2 along with the periods of two 
models from our fitting procedure, which we describe in ^3] 

3. MODEL FITTING 

We applied the pa rallel genetic-algo rithm-based fitting 
method described in Metcal fe" & Charbonneaul (120031) to 
match our white dwarf models to the observed periods of 
BPM 37093 listed in Table [T] The extension of this method 
from DBV stars to DAV stars primarily involved redefining 
the temperature range, and adding an adjustable parameter 
for the H layer mass. This led to a significant increase in 
computing time, since our code quasi-statically evolves hot 
polytropic starter models down to the temperature range of 
interest — and the DAVs are considerably cooler (~ 10,000 K) 
than the DBVs (~25,000 K). To make the problem com- 
putationally tractable, our initial exploration of DAV mod- 
els is necessarily limited compared to recent work on DBVs 
JMetcalfel2003hl). 



Table 1 

Observed and Calculated Periods for BPM 37093 



P u 

r obs 


i 


10 Mq (Pure C) 




1.10 M Q (Pure O) 




k i 


Peak 


O-C 


k 


I 


Peak 


O-C 


511.7 


28 2 


511.58 


+0.12 


27 


2 


512.50 


-0.80 


531.1 


29 2 


530.07 


+1.03 


28 


2 


530.00 


+1.10 


548.4 


30 2 


547.18 


+1.22 


29 


2 


547.32 


+1.08 


564.1 


31 2 


563.71 


+0.39 


30 


2 


564.89 


-0.79 


582.0 


32 2 


581.62 


+0.38 


17 


1 


582.27 


-0.27 


600.7 


33 2 


600.53 


+0.17 


32 


2 


600.33 


+0.37 


613.5 


19 1 


615.23 


-1.73 


18 


1 


612.61 


+0.89 


635.1 


35 2 


636.28 


-1.18 


34 


2 


636.22 


-1.12 



3.1. Defining the Parameter-Space 

The genetic algorithm explored a range of effective tem- 
peratures (Terr) between 10,000 K < T eff < 15,000 K with 
a resolution of 100 K. T his easily enco mpasses the empiri- 
cal DAV instability strip ( Bergero n et alJl2004l) . and allows 
for any possible shifts in the temperature scale of our mod- 
els due to differences in the cons titutive physics (e.g., see 
Metcalfe. Montgomery. & Kawaler 2003). Following the rec- 
ommendation of Bergeron et al. ( 1995), we fix the convective 
efficiency to the ML2/a=0.6 prescription. 

Recent evolutionary calculations suggest that the He and 
H layer masses for massive white dwarfs like BPM 37093 
should be near l ogfAfHs/M*) ~ -3.1 and log(M H /M*) ~ 
-5.8 respectively ( Altha us et alJl20 03). We allowed the ge- 
netic algorithm to fit for He layer masses between -4.0 < 
log(M He /M st ) < -2.0 with a resolution of 0.02 dex. Helium 
layers with log(Mu e /M*) > -2.0 would theoretically initiate 
nuclear burning at the base of the layer. This would likely re- 
sult in models that are too luminous to be consistent with the 
observations of BPM 37093, so we do not consider such thick 
layers. 

A similar consideration leads us to exclude H layer masses 
with log(Mn/M*) > -4.0. Our search covers the range -8.0 < 
log(M H /M*) < -4.0 with a resolution of 0.04 dex, and sub- 
ject to the constraint that the H layer be at least two orders of 
magnitude thinner than the He layer. This avoids di fficulties 
with overlapping transition zones (see lBradle\ll99l . 

Crystallization provides a new energy source for a cooling 
white dwarf: the latent heat of crystallization. This supple- 
ments the thermal energy of the ions and delays the cooling, 
so it primarily affects the ages of models at a given T e s. Al- 
though the star cools more slowly, the thermal structure at 
a given temperature is nearly identical to what it would be 
in the absence of crystallization. The change in density due 
to the transition from liquid to soli d is only a few parts per 
thousand (Lamb & Van Horn 1975), so the global mechani- 
cal structure of the star is also relatively unperturbed. While 
these evolutionary effects of crystallization have only a mi- 
nor influence on the pulsations, the presence of a solid core 
can greatly affect the pulsations. The non-radial g-modes 
are unable to penetrate the solid-liquid interface because the 
non-zero shear modulus of the solid effectively excludes their 
shear m otions, and the os cillations are confined to the fluid 
regions (Montgomery & Winget 1999). 

Since we are only interested in fitting the pulsation peri- 
ods, we can treat the degree of crystallization as an adjustable 
parameter for a given equilibrium model. We simply modify 
the inner boundary condition of the pulsational analysis for 



TESTING CRYSTALLIZATION THEORY 



3 



a given degree of crystallization. This procedure is an order 
of magnitude less computationally expensive than the evolu- 
tion part of the problem, so we can afford to repeat it several 
times. For each evolved model requested by the genetic algo- 
rithm, we calculated the pulsation periods for ten values of the 
crystallized mass fraction between 0.0 < M a /M* < 0.9 and 
returned the root-mean-square (rms) differences between the 
observed and calculated periods for the value of M cl that min- 
imized the residuals. The spherical degree (£) of the modes 
was also considered to be unknown, so we calculated all of the 
I = 1 and 1 = 2 modes between 450 and 700 s for each model. 
Since there are always more 1 = 2 modes for a given model, 
we required them to be closer to the observed period by a fac- 
tor (Ng = 2/Ni = \) to be considered a better match. In effect, we 
optimized both M a and the mode identification internally for 
each model evaluation, while the genetic algorithm optimized 
the values of the other three parameters. 

In this initial study, we performed fits for three fixed val- 
ues of the stellar mass (M t ) correspon ding to several pub- 
lished estimates. Berger on et alJ {T995) fit Balmer line pro- 
files to obtain the highest mass esti mate (M t = 1.10 Mr , 
according to a recent reanalysis by IBergeron et alJ 120041) . 
The lowest mass estimate (M* = 1.00 M ) was based 
on photometric observations and the trigonometric parallax 
(Bergeron, Leg gett. & Ruizl 120011) . An intermediate valu e 
(M» = 1.03 M ) was obtained bv iKoester & Allardl li2000). 
who used a combination of UV spectra, photometry, and the 
parallax. 

For each of these three masses we performed fits using 
models with pure C and pure O cores, which extended to 0.98 
M r /M* in fractional mass. By considering only pure core 
compositions, we do not have to consider the effect of phase 
separation of C and O during crystallization, or the subsequent 
mixin g of the remaining fluid layers (e.g., see Sala ris et"afl 
1997). Although these effects could be important, they are 
expected to have less impact for large degrees of crystalliza- 
tion, as we discuss in §0] 

3.2. Hare & Hound Exercises 

Before fitting the observed periods of BPM 37093, we per- 
formed three 'Hare & Hound' exercises to determine what 
systematic errors might arise due to our limited exploration 
of M cl , and from fixing the mass and the core composi- 
tion. We constructed a standard 1.00 M C-core model with 
T et£ = 1 1 , 700 K, log(M H e/M„) = -3.10, log(M H /M*) = -5.76, 
and M cl = 0.5 M*. We then produced three new models that 
each differed from this standard model slightly. Model 'X' 
had M a - = 0.53 M*, model 'M' had a mass of 1.03 M Q , and 
model 'C had a uniform core of 50% C and 50% O. Next, we 
calculated the periods of each of these modified models and 
tried to match them using the genetic algorithm with the mass 
fixed at 1 .00 M Q and C-core models. 

In each case, the genetic algorithm found an optimal model 
with parameters that were shifted slightly from their input val- 
ues. The optimal match for model 'X' had M a - = 0.5 M», but 
r e ff, and the He and H layer masses were offset by (-100 K, 
+0.06 dex, -0.04 dex) respectively. This suggests that our 
coarse sampling of M a - may lead to small systematic errors 
in the other three parameters. The best match for model 'C 
had the same r e ff as the input model, but the He and H layer 
masses and M CI were offset by (-0.02 dex, -0.32 dex, -0.2) 
respectively. This implies that by performing fits only with 
pure C and pure O cores, we may underestimate the actual 
value of M CI and the layer masses. The largest offsets came 



Table 2 

Fixed-Mass Optimal Models for BPM 37093 



Parameter 


1.00 Mq 


1.03 Mq 


1.10M Q 


Pure C 


Pure O 


PureC 


Pure O 


PureC PureO 


7k (K) 


13,700 


14,500 


11,100 


11,300 


10,500 11,500 


log(M H e/M»)... 


. -2.20 


-2.00 


-2.12 


-2.26 


-2.60 -2.22 


log(M H /M,)... 


-4.40 


-5.00 


-4.28 


-4.36 


-4.60 -5.64 


Mcr/M, 


0.90 


0.90 


0.90 


0.90 


0.90 0.90 


cp (s) 


1.24 


1.04 


1.08 


1.14 


0.95 0.86 



from the fit to model 'M', which overestimated the values by 
(+100 K, +0.72 dex, +0.40 dex, +0.2) for r eff , the He and 
H layer masses, and M a respectively. By fitting for only a 
few masses, we might find a spuriously high value of M CI 
and determine the layer masses only to within a factor of five. 
The above results for models 'X' and 'M' can be understood 
in terms of the calculations of Montgo mery & Wingell d 19991 
their §7.4 and Eq. 7), who showed that the average period 
spacing is much more sensitive to changes in M. t than M cr . 

3.3. Application to BPM 37093 

With these limitations firmly in mind, we applied our fitting 
method to the obse rved pulsation periods of BPM 37093. As 
mentioned in jQ.ll we repeated the fitting procedure six times: 
using pure C and pure O cores for each of three fixed masses. 
The results of these fits are shown in Table [2] In every case, 
the genetic algorithm found a model that led to rms differ- 
ences between the observed and calculated periods (ap) near 
1 s, which is comparab le to the level of accuracy we achieved 
with our DBV models ( Metcalfe 2003b). 

One of the most striking features of the model-fits in Table 
|2]is that they all have a crystallized mass fraction ofO. 9 within 
the resolution of the search (AM cr = 0.1 M*). This suggests 
that the imprint of a large value of M cr on the pulsation peri- 
ods is sufficiently strong to favor a large value of M cr in the 
fits regardless of the mass and composition. The best fits are 
achieved for the highest mass models we consider (1.10M Q ). 
The calculated periods and mode identifications for these two 
models appear with the observed periods in Tabled Note that 
an alternative identification for the pure O fit to the 582 s pe- 
riod is a (k = 31,f = 2) mode at 582.57 s, degrading the fit to 
cr P = 0.88 s. 

The values of T e s derived from spectral line profile fitting 
can provide an independent check of the asteroseismological 
mo del-f its listed in Table 13 The three mass estimates cited 
in iQ.ll also prod uced temperature estimates of: 1 1 , 550 ± 470 
K (for 1.00 Mr^ lBereeron. Leggett. & Ruizll2001l) . 11,520± 
110 K (for 1.03 : IKoester & AllardpOOOfTand 1 1, 730 ± 
200 K (for 1.10 M ? : IBergeron et alJ 1200*1 iFontaine et all 
120031) . The temperatures of the lowest mass model-fits are 
both significantly (4-6c) higher than the spectroscopic esti- 
mate, while the higher mass model-fits are mostly consistent 
with the independent measures. Our 'Hare & Hound' exer- 
cises led us to expect some small temperature shifts in the fits, 
primarily from fixing the mass and having a limited resolution 
in M cr . 

We have no independent method of determining the He and 
H layer masses, but the fit values for BPM 37093 fall between 
the canonica l values of M He ~ 10" 2 M* and M H ~ 10~ 4 M» 
( IWoodill 99% . and the theoretical values from lAlthaus et alJ 
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4. DISCUSSION & FUTURE WORK 

We have conducted the first large-scale exploration of mod- 
els to fit the observed pulsation periods of the potentially crys- 
tallized white dwarf BPM 37093. We have limited this ini- 
tial study to consider several fixed values of the stellar mass 
and the core composition to keep the problem computation- 
ally tractable. Even so, the results shown here represent more 
than 6 GHz-CPU-years of calculation, which was only prac- 
tical becaus e of the dedicated parallel computers available to 
this project (Metcal fe & Natherl2000l) . 

All of our optimal models lead to M cr = 0.9±0.1M,. Our 
'Hare & Hound' exercises suggest that we might be overes- 
timating M a - for individual fits by fixing the mass, but also 
imply that fixing the core composition could cause us to un- 
derestimate M cr by a comparable amount. The best of our six 
fits is the 1 . 10 M Q O-core model, which has an rms period dif- 
ference of only <7p = 0.86 s relative to the observed periods of 
BPM 37093. The theoretical value of M cr at the temperature, 
mass, and composition of this model is M cl (theor.) = 0.93, 
suggesting that major revisions to the existing theory may not 
be required. 

With a significant increase in computing resources, we 
should be able to extend this model-fitting method to treat 
a range of masses spanning the spectroscopic values. This 
would have the advantage of eliminating the potential sys- 
tematic errors in our fit parameters from fixin g the mass. As 
was shown by Montgomery & Winaet ( 1999), the periods of 
the modes are sensitive to very small changes in M cr , sug- 
gesting that a more thorough exploration of this parameter 
(AM C1 - ~ 0.01 M») would also be fruitful — allowing us to test 
crystallization theory with an unprecedented precision. 

Realistic stellar models do not predict pure compositions, 



but ra ther a C/O mixture that va ries as a function of radius 
(e.g., ISe gretain & ChabrierlU993l | 



alaris et aTll2000l How- 
ever, for the large degrees of crystallization which we have 
found (M cl ~ 0.9 M*), the liquid mantle above the crystallized 
core is expected to be fully mixed and significantly enriched 
in carbon due to phase separation. As a result, the original 
C/O profile in the liquid portion of the core will have been 
erased, removing any mode trapping properties which this re- 
gion would have had on the the pulsation modes. In this case, 
the mode trapping would be primarily due to the envelope 
transition zones of H and He, simplifying the problem. Thus, 
in future fits it might make sense to treat the liquid C/O mantle 
above the crystallized core as a uniform mixture of C and O, 
and to fit for the relative abundances of these two elements — 
providing additional insights about phase separation. 

In the near future we expect that the Sloan Digital Sky Sur- 
vey will u ncover several new m assive hydrogen-atmosphere 
pulsators (Mukadamet al. 2004), and each additional object 
will offer the opportunity to place independent constraints on 
crystallization theory. By sampling various masses and tem- 
peratures, we may eventually probe the equations of state of 
carbon and oxygen over a broad range of otherwise inaccessi- 
ble physical conditions. 
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